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Abstract We numerically study the evolution of magnetic fields and fluid flows 
in a thin spherical shell. We take the initial field to be a latitudinally confined, 
predominantly toroidal flux tube. For purely toroidal, untwisted flux tubes, we 
recover previously known radial-shredding instabilities, and show further that in 
the nonlinear regime these instabilities can very effectively destroy the original 
field. For twisted flux tubes, including also a poloidal component, there are 
several possibilities, including the suppression of the radial-shredding instability, 
but also a more directly induced evolution, brought about because twisted flux 
tubes in general are not equilibrium solutions of the governing equations. 
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1. Introduction 

The solar tachocline is the transition zone between the uniformly rotating ra- 
diative interior and the differentially rotating convection zone. First discovered 
helioseismically in the early 1990s, it continues to be probed to this day. Of 
particular interest is how thin it is; the angular-velocity profile changes dra- 
matically over no more than about 5% of the solar radius. This intense shear 
is one of the reasons why the tachocline is now generally believed to be the 
seat of the solar dynamo, with poloidal fields being sheared out to produce very 



strong toroidal fields. See Hughes, Rosner, and Weiss (|2007p for reviews of many 



different aspects of the tachocline, including its role in the solar dynamo. 

Given the combination of strong magnetic fields and shear, it was quickly 
realized that so-called magneto-shear instabilities could potentially play an im- 
portant role in the dynamics of the tachocline. Indeed, in other contexts, long 
before the tachocline was even discovered, it had been known that both magnetic 
fields ( |Gough and Tayler, 1966P and shear ( [Watson, 198ip separately can lead 
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to instabilities. The first to study this problem in the tachocline context were 
IGilman and Foxl (|1997p ; one very interesting result that they obtained was that 
a combination of magnetic fields and shear can be unstable even though each 
ingredient separately would be stable. IGilman and Fox, carried out a 2D calcu- 
lation, confined to the surface of a sphere, with no radial variations allowed. 
Subsequent work extended this to quasi-3D shallow-water and 3D thin-shell 
models, including also both linear onset and nonlinear equilibration studies. See 
IGilman and Cally| (|2007p for a review of this work. 

We will consider a different type of two-dimensionality, namely axisymmetric 
solutions. These have not received as much attention as some of the 3D so- 
lutions, but recent linear-onset calculations ( [Cally, Dikpati, and Gilman, 2008 



Dikpati et al., 2009) indicate that if the toroidal field is concentrated into latitu- 



dinal bands, instabilities may arise that shred the field in the radial direction. In 
this paper we consider the nonlinear evolution of such banded fields. For purely 
toroidal, untwisted flux tubes, we obtain results in good qualitative agreement 
with the linear-onset calculations. For mixed toroidal plus poloidal, twisted flux 
tubes (which were not considered before) we show that the field evolves not just 
via the onset of instabilities, but much more directly, simply because in general 
it is not an equilibrium solution of the governing equations. For twisted flux 
tubes there is then a variety of possible outcomes, depending on the strengths 
of both the toroidal and poloidal components. 



2. Equations 

The equations we wish to solve are the (Boussinesq) Navier-Stokes equation 

— + V • Vv = -Vp+ (V X a) X a+ S'e^ + eV^v, (1) 
ot 

the magnetic induction equation 

^ = V X (v X a) + eV^a, (2) 

and the entropy equation 

— +vWS^ -N^v ■ er + eV^S. (3) 

Here v and a denote the fluid and Alfven velocities, respectively, and S the 
entropy. N is the Brunt- Vaisala frequency, assumed constant. Length has been 
scaled by the inner edge of the tachocline, so we will be interested in solving 
this system in the interval r g [1, 1.05]. Time has been scaled by the equatorial 
rotation frequency, v and a as length/time. 

Except for the inclusion of the diffusive terms (eV^-), these equations are the 
same as in |Gally, Dikpati, and Gilman| (|2008p . hereafter referred to as CDG08. 
Some diffusivity must be included here to ensure numerical stability, but values 
in the range e = 10~^ to 10~^ yielded similar results, indicating that diffusivity 
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is not significantly affecting the evolution. The results presented here are all at 
e = 2 X 10"^. 

While the equations may be much the same, the subsequent analysis is very 
different from that of CDG08. They considered only the linear onset of instabil- 
ity, and only for high radial wavenumber (fc), thereby ultimately eliminating r 
entirely, instead simply having as a parameter in the equations. In contrast, we 
are interested in a direct numerical solution, in the finite interval r £ [1, 1.05], 
allowing an arbitrary radial dependence, and including also the full nonlinear 
evolution of the solutions. 

For axisymmetric solutions, it is convenient to decompose v and a as 

V = -ye^ + V X (V'e^), a = 6e<j + V x (ae^). (4) 

The toroidal parts, v and 6, are the azimuthal components of the given vec- 
tors; the poloidal parts, V and a, are the streamfunctions of the meridional 
components. The original Equations (1) and (2) then become 

^^v = Pi{b,a)-Pi{v,^)+eD''v, (5) 

g-^D^^ = -9gS + P2{b, b) + P2p'a, a) - P^iv, v) - P2{D^, ^) + eD^, (6) 
^b = P2{v,a)-P2{b,^)+tD\ (7) 



where 



and 



^a = Pi(V',a)+e£>V (8) 



= V" - l/(rsine)", (9) 



Pi(X,y) =e^- [(V X (Xe^)) X (Vx (Fe^))], (10) 

P2(X, F) = • V x [{Xe^) x (V x (Ye^))] . (11) 
We then wish to solve (5) and (6), with stress- free boundary conditions 

§-,{^=^ = ^2^ = ^ at r = 1,1.05, (12) 
(7) and (8), with perfectly conducting boundary conditions 
d 

— (br) =0 = at r = 1, 1.05, (13) 

or ^ ' 

and (3), with 5 = at r = 1, 1.05. These boundary conditions are somewhat 
artificial, but any boundary conditions would necessarily be artificial. To prop- 
erly capture all of the dynamics of the tachocline, it should not be studied in 
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isolation, but rather as part of a global model that also includes the interior 
and the convection zone. However, while there are models that aim in this 
direction ( |Brun and Toomre, 2002^ , they are so complicated that one cannot 
focus specifically on aspects such as tachocline instabilities. To study tachocline 
instabilities, one must therefore adopt simplified models such as ours, despite 
the inevitably unnatural boundaries where the real Sun has none. 

We solve these equations using the numerical code described bv IHoUerbachl 
POOOP . in which the radial structure is expanded in terms of Chebyshev poly- 
nomials, the angular structure in terms of spherical harmonics, and the time- 
stepping is done via a second-order Runge-Kutta method. Resolutions ranging 
from 50 x 1600 to 80 x 2400 in (r, 9) were used, and were all checked to ensure 
that the results were adequately resolved. A timestep of 10"'^ was sufficiently 
small to ensure stability. 

To facilitate comparison with CDG08, we choose initial conditions much the 
same as theirs. The initial entropy is simply S* = 0, so any buoyancy effects 
arise entirely out of the subsequent evolution. For the flow, we take the solar- 
like differential rotation profile uj = v/irsvuO) = 1 — O.18cos^0, where we note 
that (1) is written in a non-rotating frame, so w here must include an overall 
rotation, not just the Pole-to-Equator differential rotation — 0.18cos^ 9. 

For the toroidal field, we take 

b = Ap[e-4(^-'^)'/^'(i-'^') - e-4(M+d)Vv^^(i-d^)] sin^/,., (M) 

where /i — cos 9. That is, b consists of two oppositely directed bands in the two 
hemispheres, with position d, latitudinal bandwidth W ^ and amplitude A. The 
normalization factor p is adjusted such that the maximum field strength max£/(6) 
is A/2. 

Values of A up to one, corresponding in dimensional terms to max5i(&) « 10^ 
G, are of greatest interest in the solar context, but other (especially younger) 
solar-type stars may well have even larger values. We will present runs for A — 1 
and 2. A range of possibilities for d and W were considered, and yielded qualita- 
tively similar results. We therefore show results only for d = 0.5, corresponding 
to bands situated ±30° from the Equator, and W = 7r/36, corresponding to a 
bandwidth of 5°. 

Thus far our initial conditions are exactly as in CDG08; see also Dikpati and GilmanJ 



(|1999p . who were the first to introduce banded toroidal fields of this type. To 
this toroidal field, we now add the poloidal field 

a = AV[e-4(^-'^)'/^'(i-'^') + e-4(M+d)Vv^^(i-rf^)] sin^(r - l)(r - 1.05). (15) 

The poloidal field thus has the same banded structure as the toroidal, but a is 
equatorially symmetric, whereas b is anti-symmetric. Both of these symmetries 
correspond to the standard "dipole" solutions of solar dynamo theory. The dif- 
fering radial dependencies, (r — l)(r — 1.05) for a versus 1/r for 6, are dictated 
by the different boundary conditions (13) that a and b are supposed to satisfy. 
The normalization factor p' is adjusted such that maxe(a) is A' . 

The amplitude A' = fA, so the factor / gives the ratio of poloidal to toroidal 
fields. In addition to the f — 0, untwisted flux tubes, we will take / = ±10"* 
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Figure 1. Prom left to right, contours of b, uj, ip, and S, with contour intervals of 0.2, 0.02, 
5 X 10~^, and 0.1, respectively. In the plots for tf> and S, grey denotes positive values, white 
negative (■(/>> corresponds to clockwise circulation). From top to bottom, the three rows are 
for = 10, 10^, 10^. A = 1 and time t = 5 for all three. The range of r is 1 to 1.05, but has 
been stretched by a factor of 10, and hence looks like 1 to 1.5. 

and ±4 x 10^^ for the twisted flux tubes. The sign of / determines whether 
the tubes are twisted in a left- or right-handed sense. It is not certain which 
is more appropriate to the tachocline ( [Fan, 2004| , so we consider both possi- 
bilities, and show that they yield qualitatively similar behavior. Regarding the 
amplitudes of /, for A = 1 and |/| = 10~^, the maximum values of {ar,ag,a^) 
are (0.0022, 0.0092, 0.5). Even for |/| = 4 x IQ-* the field is thus predominantly 
azimuthal, with only around one twist over the full circumference </> £ [0,27r]. 
We will see though that even such relatively weak poloidal fields as this can 
significantly influence the subsequent evolution. 

3. Results 

3.1. Untwisted Flux Tubes 

We begin by considering the influence of the stratification. Figure 1 shows results 
for = 10, 10^ 10^, all for A = 1 and / = 0. Within each row, the first panel 
shows the toroidal field b, the second panel the angular velocity uj = v/{rsmd), 
the third the meridional circulation ip, and the fourth the entropy S. (According 
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Figure 2. Prom left to right, contours of 6, oj, ip, and 5, with contour intervals of 0.2, 0.05, 
5 X 10~®, and 1, respectively. In the plots for ui, grey indicates regions where a; > 1; for ^ 
and 5, grey denotes positive values, as in Figure 1. AT^ = 10^, A = 2, and from top to bottom 
t = 2.5, 5, and 10. 

to Equation (8), if a = initially, it will remain zero.) Only the range 9 G 
[45°, 75°] is shown, centered on the flux tube at = 60°. The radial direction 
has been stretched by a factor of 10, that is, the actual gap r € [1, 1.05] has been 
stretched to look like [1, 1.5]. 

Turning to the variation with N'^, we see that for = 10 the solution has 
changed significantly from its initial condition, whereas for A^^ = 10^ it is almost 
unchanged. The reason for this is easy to understand: If initially only b and v are 
non-zero and if a is always zc;ro then according to Equations (5) and (7), the 
only way (apart from the very weak dissipation) for either 6 or z; to change is by 
inducing a meridional circulation ip. Now, according to Equation (6), the terms 
P2(&, b) — P2(v, v) will indeed induce such a circulation; these terms are zero only 
if -^[b"^ — v^) = 0, which in general is not the case. However, for increasingly 
strong stratification, the buoyancy term (?'~^ff ) can very effectively suppress 
the tendency to drive a meridional circulation. That is, having only b and v non- 
zero is not quite an equilibrium solution to the governing equations, but if 
is sufficiently large, only a very weak circulation ip will be induced, so according 
to Equations (5) and (7), b and v will remain almost unchanged. This justifies 
linear stability analyses such as those of CDG08, where a basic state is simply 
imposed for b and v. In the remainder of this paper we will consider only the 
strongly stratified case N'^ = 10^. 
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Figure 3. From left to right, contours of b, a, u), ip, and S, with contour intervals of 0.2, 10~^, 
0.1, 2 X 10^"', and 0.2, respectively. In the plots for w, grey indicates regions where w > 1; for 
tl) and S, grey denotes positive values, as in Figure 1. = 10^, A = 1, and / = 4 X lO"", 
corresponding to clockwise circulation for the poloidal field. The top row is at t = 2.5, the 
bottom row at t = 5. 
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Figure 4. Prom left to right, contours of 6, a, lu, i/>, and S, with contour intervals of 0.2, 10 
0.1, 2 X 10-^, and 0.2, respectively. = 10^, A = 1, and / = -4 x 10"*, corresponding to 
counter-clockwise circulation for the poloidal field. (Incidentally, it is perhaps explicitly worth 
noting that because of the ten-fold stretching in the radial direction, ag is actually ten times 
greater than the spacing of the contour levels might suggest. So in fact ag is greater than Or.) 
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Figure 5. Prom left to right, contours of 6, a, u>, tj}, and 5, with contour intervals of 0.2, 
5 X 10~^, 0.1, 5 X 10~^, and 1, respectively. In the plots for oj, grey indicates regions where 
u) > 1; for ip and S, grey denotes positive values. AT^ = 10^, A = 2, and / = 10"''. Prom top 
to bottom t = 2.5, 5, and 10. 



Figure 2 shows the effect of doubhng the field strength, to A ~ 2. Now even 
the strong stratification is not enough to stabihze the solution. Instead, we see 
the development of precisely the radial-shredding instabilities previously studied 
by CDG08. Here though we follow the full nonlinear evolution, and discover that 
by t = 10 the instability has almost completely obliterated the original flux tube. 

3.2. Twisted Flux Tubes 

Figures 3 and 4 show the results for A = 1 and / = ±4 x 10"**. That is, the 
amplitude A = 1 is as in Figure 1, too low for the shredding instability to occur, 
and indeed it doesn't. Nevertheless, the solutions also do not remain virtually 
stationary, as in the bottom row in Figure 1. Instead, we see the development 
of highly localized jets in the angular velocity u). To understand their origin, we 
return to Equation (5), where the term Pi (6, a) is clearly responsible; this term 
is zero only if contours of b and a coincide, which in general is not the case. 
Furthermore, because there is no buoyancy force in this equation, no amount of 
stratification can suppress this effect, very much unlike Figure 1. 

Comparing Figures 3 and 4 in detail, we note that reversing the sign of a (that 
is, the sense of twist in the tube) has exactly the effect one might expect; revers- 
ing the sign of Pi (6, a) simply reverses the jets. Otherwise the evolution is much 
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Figure 6. Prom left to right, contours of 6, a, w, V", and S, with contour intervals of 0.2, 
5 X 10~^, 0.1, 5 X 10~^, and 1, respectively. In the plots for oj, grey indicates regions where 
u) > 1; for ip and S, grey denotes positive values. AT^ = 10^, A = 2, and f = Ax 10~*. Prom 
top to bottom t = 2.5, 5, and 10. 



the same, and in both cases the flux tubes largely maintain their strength. Note 
also that jets as concentrated as this would not be detectable by helioseismology, 
so phenomena such as these could conceivably exist in the real tachocline. 

Figure 5 shows results for A = 2 and / = lO"**. The toroidal field is therefore 
as in Figure 2, whereas the poloidal field is half as strong as in Figure 3 (so the 
nonlinear term Pi (6, a) is just as strong as in Figure 3). Initially we see the same 
jets as in Figures 3 and 4 (and reversing the sign of / again merely reverses the 
jets). By t = 5 the evolution is dominated by the same shredding instability as 
in Figure 2, and the final result is much the same, with the original fiux tube 
largely destroyed. Poloidal fields as weak as this therefore have relatively little 
influence. However, if we increase the poloidal field strength to / = 4 x 10~*, 
it does have a very significant influence, as illustrated in Figure 6. The radial- 
shredding instability is now completely suppressed, and the flux tube persists up 
to t = 10 (and beyond). Note also how the solution has adjusted itself so that 
contours of h and a do now largely coincide, and correspondingly these localized 
jets are greatly reduced in strength. 
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4. Conclusions 

We have considered the evolution of flux tubes in a thin spherical shell, intended 
to model solar-type tachoclines. For untwisted tubes we obtain the same radial- 
shredding instabilities that have previously been studied in the linear regime, and 
show that in the nonlinear regime these instabilities can very efficiently destroy 
the original flux tube, simply by shredding it to sufficiently short length-scales 
for it to dissipate. 

For twisted flux tubes, there are a number of possibilities. If the toroidal fleld 
is too weak for instabilities to set in, the solution will nevertheless evolve, via 
the formation of differential-rotation jets, driven directly by the Lorentz forces 
associated with the twist in the tube. These jets induce a certain amount of 
structure in the flux tubes, but not enough to disrupt it as the instabilities did. 

If the toroidal field is sufficiently strong, and the poloidal field very weak, 
the shredding instabilities develop much as before, and simply overwhelm the 
jets driven by the twist. Finally, if the poloidal field is somewhat stronger - but 
still much weaker than the toroidal - it can suppress the shredding instability. 
Twisted flux tubes can therefore exist at considerably greater fleld strengths 
than untwisted tubes. 

Future work will extend this model to 3D and study the interaction of some of 
the effects presented here with some of the previously known non-axisymmetric 
magneto-shear instabilities ( Gilman and Cally, 2007[ ). 
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